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I.  STATEMENT  OF  THE  PROBLEM  STUDIED 

Flows  over  helicopter  blades  are  very  complex  turbulent  flows.  The  blades  experience 
dramatically  different  flow  field  at  various  azimuthal  angles.  Also,  a  rotor  normally  consists  of 
many  elastic  blades,  with  a  strong  coupling  between  aerodynamics  and  structure.  The  problem  is 
indeed  multidisciplinary.  Current  helicopter  blade  designers  use  computational  models,  which 
depend  heavily  on  experimental  data  and  cannot  be  used  to  predict  any  novel  design,  which  is  a 
significant  departure  from  existing  designs.  To  simulate  this  multiscale  and  multidisciplinary 
physics  with  confidence,  we  have  developed  a  robust  multidisciplinary  computational  tool  called 
WINDUS-HELI  based  on  the  WIND-US  code  by  coupling  CFD  with  CSD,  adopting  state-of- 
the-art  numerical  approaches,  and  applying  high  fidelity  physics  models.  Very  reasonable  results 
have  been  obtained  for  both  the  aerodynamic  loads/performance  and  acoustics  predictions  for  all 
the  validation  cases  studied. 

II.  SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 
WINDUS  Enhancement 

The  WIND-US  code  was  chosen  for  this  study  due  to  the  fact  that  it  has  the  capability  of  solving 
the  Navier-Stokes  equations  on  a  computational  domain  that  is  discretized  using  an  arbitrary 
combination  of  structured  and  unstructured  grids.  [Mani  et.  al.]  In  addition,  the  code  can  handle 
possible  overlapping  grids  in  the  interior  of  the  domain,  which  is  very  desirable  for  hybrid 
URANS/LES  simulations  and  for  the  development  of  novel  blade  shapes.  In  this  work,  the 
enhancements  made  for  WINDUS-HELI  code  include  second-order  temporal  accuracy  and 
global-Newton-iterations  with  dual-time-stepping,  moving-grid  capability  for  rotation  and 
translation  of  rotor  [Jain  et.  al.],  unstructured-mesh-deformation  for  controls  motion  and  elastic 
deformation  of  the  rotor  blades,  parallel  efficiency  improvement,  multi-grid  scheme,  anisotropic 
non-linear  k-s  model  [Shih  et.  al.],  DES  model,  PRNS  models  [Shih  et.  al.]  and  a  free  wake 
model.  Interfaces  have  also  been  developed  in  WINDUS-HELI  for  RCAS  [Saberi  et.  al.]  and 
PSU-WOPWOP  codes  for  flow  structure  interaction  and  acoustics  signature  predictions.  [Figure 
la]  A  sleep-awake  mechanism  was  implemented  in  WINDUS-HELI  for  loose  coupling  with 
RCAS.  [Figure  lb]  The  loose  coupling  process  starts  from  RCAS  that  trims  the  model  with 
internal  aerodynamics.  After  multiple  periodic  steady  state  calculations,  RCAS  outputs  the  first 
blade  motion  file  and  creates  the  flag  to  inform  WINDUS-HELI  that  the  motion  file  is  ready  and 
goes  to  sleep  waiting  for  the  airloads  flag.  The  WINDUS-HELI  master  processor  reads  the 
motion  file,  updates  the  cubic  spline  coefficients  used  for  interpolation  and  broadcasts  this 
information  to  all  the  slave  processes  and  starts  the  CFD  calculation.  For  a  4-bladed  rotor,  after 
finishing  a  quarter  revolution,  WINDUS-HELI  integrates  the  airloads  and  outputs  the  airloads  at 
the  specified  radial  locations,  creates  the  flag  to  tell  RCAS  that  the  airloads  file  is  ready  and  goes 
to  sleep  waiting  for  the  motion  flag.  When  RCAS  detects  the  flag,  it  wakes  up,  deletes  the  flag, 
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Figure  1  a)  WINDUS-HELI  code  Enhancement,  b)  automated  data  exchange  mechanism 

reads  the  updated  airloads.  RCAS  calculates  the  difference  between  the  CFD  airloads  and  RCAS 
airloads  from  the  previous  trim,  and  assigns  the  difference  to  the  mechanical  loads.  RCAS  re¬ 
trims  with  the  delta  force  and  outputs  the  updated  motion  data  and  creates  the  motion  flag. 
WINDUS-HELI  detects  the  flag,  wakes  up,  deletes  the  flag,  reads  the  motion  data  and  starts  the 
CFD  calculation.  The  whole  process  is  repeated  until  there  is  little  change  in  the  RCAS  and 
WINDUS-HELI  solutions,  and  convergence  is  observed  in  the  RCAS  trim  variables,  which 
consist  of  collective  angle,  lateral  cyclic  angle  and  longitudinal  cyclic  angle. 


Aeroload  Validation 


The  validation  of  normal  force  contour  plots  of  UH-60A  cases  including  C8534,  C9017,  C9310, 
C9121,  and  C9505  are  shown  in  Figures  2-6. 


270 


(a)  (b) 

Figure  2  Contour  plots  of  normal  force  for  C8534:  (a)  Flight  test,  (b)  WIND-US  RCAS  coupling. 
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(a)  (b) 

Figure  3  Contour  plots  of  normal  force  for  C9017:  (a)  Flight  test,  (b)  WIND-US  RCAS  coupling. 
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(a)  (b) 

Figure  5  Contour  plots  of  normal  force  for  C9121:  (a)  Flight  test,  (b)  WIND-US  RCAS  coupling. 
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Figure  6  Contour  plots  of  normal  force  for  C9505:  (a)  Flight  test,  (b)  WIND-US  RCAS  coupling. 
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The  airload  comparison  for  C8534  at  radial  locations  are  shown  in  Fig.  7.  All  the  results  show 


(a)  (b)  (c) 

Figure  7  Airload  comparison  for  C8534  at  55%,  86.5%  and  96.5%  of  radius,  (a)  normal  force, 

(b)  pitching  moment,  and  (c)  chord  force. 


Acoustics  Validation 

An  acoustics  interface  has  been  developed  in  WINDUS-HELI  for  performing  acoustics  signature 
calculations.  The  interface  has  been  developed  for  the  acoustics  code,  PSU-WOPWOP  [Hennes 
et.  al.].  Two  options  are  available  for  acoustics  computations:  (a)  Impermeable  Surface  (b) 
Permeable  Surface.  User  defines  provide  a  surface  on  which  flow  quantities  are  interpolated  and 
written  out  as  a  sequence  of  time  series.  For  impermeable  surface  option,  loading  (pressure)  data 
is  computed  on  the  blade  surface  as  shown  in  Figure  8.  For  permeable  surface  option,  flow-field 
data  (density,  momentum,  and  energy)  are  computed  on  a  surface  away  from  the  blade  surface  as 
shown  in  Figure  9. 
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Figure  8  Impermeable  surface  approach  for  acoustics  calculations.  Time  series  loading/pressure 

data  is  computed  on  the  rotor  blade  surface. 


Figure  9  Permeable  surface  approach  for  acoustics  calculations.  Time  series  flow-field  data  is 
computed  on  a  surface  away  from  the  rotor  blade. 

It  was  found  that  the  permeable  surface  results  were  very  sensitive  to  the  interpolation  error  and 
therefore  special  procedure  was  developed  to  minimize  these  errors.  For  the  robustness  of  the 
interpolation  procedure  for  the  permeable  surface  option,  the  relative  motion  between  CFD  mesh 
and  acoustics  surface  should  be  minimized.  Abrupt  variations  in  CFD  mesh  size  (non-smoothed 
CFD  mesh)  can  lead  to  large  interpolation  errors.  It  was  found  that  such  interpolation  errors  in 
transferring  CFD  solution  to  acoustics  surfaces  can  lead  to  spurious  spikes  in  acoustics  pressure 
prediction  Therefore,  we  have  implemented  lead-lag,  flapping,  and  pitching  motion  of  the 
acoustics  surface  which  causes  the  permeable  surface  to  follow  the  CFD  mesh  very  closely  and 
this  minimize  the  relative  motion  between  the  two.  Based  on  the  blade  tip  motion,  collective, 
1/rev  and  2/rev  motion  are  computed  and  supplied  as  input  to  WINDUS-HELI.  The  permeable 
surface  is  positioned  according  to  these  motion  inputs  every  time  step.  The  acoustics  surface 
remains  rigid  (no  flexing,  only  rigid  body  motions).  It  was  found  that  the  artificial  were  not 
encountered  if  the  permeable  surface  is  moved  as  discussed  above. 

Acoustics  results  for  the  case  DNW1310  case  have  been  validated.  The  acoustics  calculations 
were  performed  at  the  nine  microphone  locations  using  both  the  permeable  and  impermeable 
options  and  compared  with  the  wind-tunnel  measurements  supplied  by  Sikorky.  The  permeable 
surface  is  shown  in  Figure  9.  The  comparison  is  shown  in  Figure  10.  Both  the  permeable  and 
impermeable  surface  results  predict  the  negative  peak  as  well  as  the  phase  very  well  at  all  the 
location.  The  permeable  surface  option  leads  to  better  prediction  of  the  negative  peak  at 
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microphone  location  #7.  This  could  be  due  to  the  presence  of  greater  high-speed  impulsive  noise 
component  at  this  location  which  is  known  to  be  better  captured  by  the  permeable  surface  option. 
The  positive  peak  is  not  well  captured  by  either  option.  The  reason  for  this  behavior  is  not 
understood  at  this  time.  The  effect  of  spatial  and  temporal  resolution  of  the  CFD  data  as  well 
comparison  of  the  rotating  permeable  surface  (where  the  surface  rotates,  flaps,  lead-lags,  and 
pitches  along  with  the  blade)  versus  non-rotating  surface  (where  the  surface  encompasses  the 
entire  rotor  and  does  not  move  with  the  rotor)  needs  to  be  examined  to  see  their  affect  on 
predicting  the  positive  peak.  Overall  very  satisfactory  agreement  with  the  experimental  results 
was  found.  The  permeable  surface  gives  slightly  better  prediction  in  general  for  this  case. 


Atouiiics  PtttUttKHi:  DNW 1 110,  ytiL'Ji  i 


ArtMOik*  fTcdkilda:  &NWt3IO, 


AttMdie*  PnftlUMttai:  DNW 1 110,  Mkfl 


Figure  10  Permeable  Surface  (red)  and  Impermeable  Surface  (blue)  acoustics  prediction  for 
case  DNW1310  at  microphone  location  1  thru  9  (in  order).  Experimental  measurement  is 

represented  by  black  curve. 

Acoustics  validation  of  “blind”  test  cases 

In  addition  to  the  DNW1310  cases,  acoustics  calculations  were  made  for  three  UH60A  level 
flight  cases  list  below 

(1)  C9310 

(2)  C9121 

(3)  C9505 

Acoustics  data  was  reported  in  the  form  of  time  history  data  and  frequency  data.  The  time  history 
data  (total  acoustics  pressure)  was  reported  in  Pascals  for  a  full  rotation  at  1024  locations.  The 
frequency  data  was  reported  as  a  single  sided  power  spectrum  in  dB’s  relative  to  20  pPa.  The 
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power  spectrum  is  defined  as  the  power  spectral  density  times  the  frequency  bin-width  of  the 
Fourier  analysis.  The  data  was  reported  on  a  hemisphere  of  radius  100  ft  centered  at  the  rotor  hub 
(see  Figure  11)  every  5  degrees  in  theta  and  phi.  The  theta  angle  defines  the  angle  under  the 
vehicle  and  varies  from  5  deg  (nose)  to  175  deg  (tail)  in  5  degree  increments.  The  phi  angle 
defines  the  sideline  angle.  It  varies  from  -90  deg  (starboard)  to  90  deg  (port)  in  5  degree 
increments. 


Figure  11  Hemisphere  of  radius  100  ft.  below  the  rotor  centered  at  the  rotor  hub  for  computing 

acoustics  data  for  the  "blind"  test  cases. 

The  data  was  reported  to  DARPA  for  an  independent  evaluation  of  the  acoustics  prediction 
against  the  flight  test  data.  Acoustics  prediction  for  the  three  cases  at  a  point  three  rotor  radius 
distance  from  the  hub  directly  ahead  is  shown  in  Figure  12. 


Figure  12  Acoustics  pressure  at  a  test  point  directly  ahead  of  the  rotor  in  the  rotor  plane  at  a 
distance  of  three  rotor  radius.  Cases  C9310,  C9505,  and  C9121  (from  left  to  right  above). 


Turbulence  Modeling 

The  focus  of  this  effort  was  to  introduce  state  of  the  art  turbulence  modeling  into  the  WINDUS- 
HELI  code.  When  the  program  began  the  only  appropriate  turbulence  model  available  in 
WIND-US  was  the  Spalart-Allmaras  model.  By  the  end  of  the  program  the  code  had  evolved  to 
a  point  where  the  user  can  chose  between  various  classes  of  detached  eddy  simulation  models 
(DES),  a  non-linear  k-e  model  (NLKE),  and  various  forms  of  large  eddy  simulation  models 
(LES).  The  various  models  were  validated  in  a  variety  of  flows  including  flow  over  a  circular 
cylinder,  flow  in  a  diffuser,  static  airfoil  stall,  dynamic  airfoil  stall,  hovering  rotor,  rotor  in 
forward  flight,  rotor  in  forward  flight  including  dynamic  stall.  The  turbulence  models  available 
in  WINDUS-HELI  are  as  follows: 

a)  Spalart-Allmaras 

b)  Detached  eddy  simulation  (DES) 

c)  Partially-resolved  Navier-Stokes  (PRNS) 

d)  Detached  partially-resolved  Navier-Stokes  (DETPRNS) 

e)  Non-liner  k-e  model  (NLKE) 
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f)  Smagorinsky  subgrid-scale  (SGS)  model  for  large  eddy  simulation  (LES) 

g)  Dynamic  Smagorinsky  (SGS)  model  for  LES  (DLES) 

f)  Capability  for  assigning  distinct  models  to  individual  zones,  most  often  used  with  zonal  LES 
(ZLES) 

Flow  over  a  circular  cylinder  was  used  as  an  initial  validation  for  the  various  turbulence  models. 
The  topology  of  the  wake  is  an  indication  of  the  range  of  turbulent  motions  that  can  be  simulated 
directly  with  the  hybrid  (DES,  PRNS)  and  LES  models  and  the  Strouhal  number  provides  a 
quantitative  measure  of  the  model  performance.  Figure  13  shows  a  comparison  of  the  cylinder 
wake  flows  for  the  Spalart-Allmaras  and  DES  models.  The  increased  range  of  resolved  scales  in 
the  case  of  the  DES  model  is  evident. 


Spalart-Allmaras  (SA)  URANS 


Detached  Eddy  Simulation  (DES) 


Figure  13  Comparison  of  a  circular  cylinder  wake  flow  using  the  Spalart-Allmaras  and  DES 

models. 

The  next  validation  step  was  static  and  dynamic  airfoil  stall.  The  experimental  data  of 
McCroskey  et  al.  were  used  as  a  basis  for  comparison.  Figure  14  shows  a  set  of  static  stall 
results  for  a  NACA  0012  airfoil  at  Mach  0.3.  All  turbulence  models  tested  were  able  to 
reproduce  the  stall  with  good  accuracy.  The  dynamic  stall  cases  proved  to  be  much  more 
difficult  and  accurate  results  were  obtained  only  after  systematic  studies  of  the  effects  of  wind 
tunnel  walls,  time  step,  iteration  tolerance,  numerical  dissipation  parameters,  and  turbulence 
model.  Figure  15  shows  a  summary  of  the  final  results  for  a  deep  dynamic  stall  case  with 
alpha=15+10sin(wt)  and  reduced  frequency  k=0.1.  The  Spalart-Allmaras  and  DES  models 
proved  to  be  reasonably  accurate,  whereas  the  PRNS  model  proved  to  be  less  accurate.  These 
results  were  obtained  with  a  mesh  that  faithfully  represented  the  wind  tunnel  walls.  The  time 
step  was  set  to  1 /3600th  of  an  oscillation  period. 
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NACA0012,  static,  M  =  0.3 


Figure  14  Static  stall  of  a  NACA  0112  airfoil  at  M=0.3. 


2.5 

2 

1.5 

1 
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a 

Figure  15  Dynamic  stall  of  a  NACA  0112  airfoil  at  M=0.3,  alpha=15+10sin(wt)  and  reduced 

frequency  k=0.10. 


20  Newton  iterations  were  used  with  10  subiterations  .  The  smoothing  (numerical  dissipation) 
parameters  were  set  at  (0.1, 1.0, 5.0).  The  solutions  were  found  to  degrade  if  any  of  the 
parameters  quoted  above  were  adjusted  in  a  direction  that  would  increase  numerical  error  (such 
as  increasing  the  time  step  or  increasing  the  numerical  dissipation). 

The  various  turbulence  models  were  also  tested  in  the  hovering  rotor  case  of  Caradonna  and 
Tung.  Figure  16  shows  a  typical  comparison  for  the  blade  loading  at  the  80%  radial  station.  The 
results  were  found  to  depend  very  little  on  the  turbulence  model. 
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Figure  16  Pressure  distribution  on  the  Caradonna  and  Tung  hovering  Rotor  at  the  80%  radial 

station. 


The  non-linear  k-e  model  was  developed  late  in  the  program  and  was  not  subjected  to  the 
validation  tests  discussed  above.  It  was  found  to  give  more  accurate  results  for  the  UH-60  test 
cases  however.  By  virtue  of  the  non-linear  terms,  this  model  is  better  able  to  account  for  the 
effects  of  swirl,  which  enables  a  much  more  accurate  capturing  of  the  blade  tip  vortices.  This 
effect  is  shown  in  Figure  17  where  the  computed  pressure  distribution  near  the  wing  tip  is 
compared  with  the  measurements  of  Chow,  Zilliac,  and  Bradshaw.  This  plot  shows  that  either 
the  non-linear  k-e  model  or  the  Spalart-Allmaras  model  with  a  rotational  correction  is  required  to 
obtain  the  correct  solution  in  the  vicinity  of  the  vortex.  Visualizations  of  the  vortex  itself  show  it 
to  be  resolved  best  when  the  nonlinear  k-e  model  is  used.  Thus  this  model  was  adopted  for  all  of 
the  UH-60  test  cases. 


z/b-0.8889 


Figure  1 7  Pressure  distribution  very  near  the  tip  of  the  Chow,  Zilliac,  and  Bradshaw  wing.  The 
non-linear  k-e  model  is  referred  to  as  the  cubic  k-e  model  in  the  figure  legend. 
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Vortex  Wake  Model 

We  have  also  implemented  a  conventional  vortex  model  where  the  Navier-Stokes  computation  is 
limited  to  within  2-3  chords  of  from  the  blade  surface  and  the  vortex  wake  model  alone  is  used 
for  the  far  field.  The  sketch  of  the  arrangement  is  shown  in  Figure  18.  The  vortex  circulation  is 
determined  as  a  function  of  time  from  the  computed  blade 


Figure  18  Schematic  of  the  vortex  wake  model. 

loads.  The  vortex  wake  is  then  discretized  into  a  large  number  of  linear  segments,  which  move 
in  time  in  such  a  way  that  the  system  is  always  force-free.  The  velocities  induced  by  the  vortex 
system  are  fed  into  the  boundary  conditions  applied  on  the  outer  surface  of  the  WINDUS-HELI 
zone.  Figure  19  shows  the  Navier-Stokes  domain  for  a  two-bladed  NACA-0012  hover  test  case. 
The  outer  boundary  is  placed  at  a  distance  of  two  chords  away  from  the  blade  surface.  The 
computed  solutions  are  compared  with  the  hover  measurements  of  Caradonna  and  Tung  .  The 
blades  are  set  at  a  collective  pitch  angle  of  8  degrees  and  the  tip  Mach  number  is  0.52. 


Figure  19  Navier-Stokes  domain  for  the  two-bladed  test  case. 


Figure  20  shows  a  comparison  of  the  pressure  computed  pressure  distributions  at  several 
different  radial  stations.  Results  are  shown  for  a  “full  domain”  case  where  the  vortex 
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Figure  20  Comparison  of pressure  distributions  for  the  NACA-0012  hover  case  using  the  vortex 
wake  model.  From  top  to  bottom  the  figures  and  left  to  right  are  for  r/R=0.50,  0.68,  0.80,  0.89, 

and  0.96. 


wake  is  computed  directly,  a  case  on  the  small  or  “partial  domain”  without  a  vortex  model  and 
finally  a  case  on  the  partial  domain  with  a  vortex  model.  As  expected,  the  results  for  the  case  of 
the  partial  domain  without  a  wake  model  (green  curves)  show  poor  agreement  with  the 
experimental  data.  The  results  are  improved  greatly  for  all  stations  when  the  wake  model  is 
added  (blue  curve).  Finally,  the  results  with  the  wake  model  also  agree  well  with  the  simulation 
done  on  the  full  domain  where  the  wake  vortices  are  computed  directly. 

High-order  Spectral  Difference  approach 

In  order  to  tackle  the  following  challenges  in  the  helicopter  flows — interactions  between  blades, 
vortex  and  wakes,  complex  and  dynamic  geometries,  and  disparate  length  scales — a  high  order 
multi-domain  spectral  difference  method  has  been  developed  for  the  three  dimensional  Navier- 
Stokes  equations  on  unstructured  hexahedral  grids.  The  code  has  been  verified  with  the  property 
of  high  order  accuracy  and  spectral  convergence  in  space;  in  time  integration,  the  code  includes 
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both  explicit  and  implicit  schemes.  An  efficient  implicit  lower-upper  symmetric  Gauss-Seidel 
(LU-SGS)  solution  algorithm  has  been  developed.  The  LU-SGS  solver  is  preconditioned  by  the 
block  element  matrix,  and  the  system  of  equations  is  then  solved  with  a  LU  decomposition 
approach.  The  implicit  solver  has  shown  more  than  an  order  of  magnitude  of  speed-up  relative  to 
the  multi-stage  Runge-Kutta  time  integration  scheme  for  several  demonstration  problems.  The 
high  order  curve  boundary  handling  is  also  included  in  the  code,  and  the  huge  difference  between 
low  and  high  order  boundary  is  demonstrated  by  simulating  an  inviscid  flow  over  a  sphere.  The 
code  has  been  parallelized  and  successfully  used  to  run  the  vortex  ring  case  and  the  hovering 
rotor  case,  which  also  demonstrates  that  the  code  can  handle  curved  boundary  with  higher  order 
curve  surface  fitting,  possess  high  resolution  of  vortices,  run  the  moving  grid,  and  deal  with 
complex  geometries.  A  NACA  0012  two  bladed  rotor  [Caradonna,  F.X.,  and  Tung,  C.]  has  been 
computed  with  SD  code.  More  than  three  revolutions  of  tip  vortex  trajectories  have  been 
observed  and  the  potential  of  wake  capturing  has  been  demonstrated.  [Figure  21]  The  spectral 
difference  code  has  been  integrated  into  the  WINDUS-HELI  code,  which  has  been  tested  on 
some  simple  problems. 


Figure  21  Iso-surface  of  down  wash 


Summary 

In  summary,  the  WIND-US  based  helicopter  rotor  code  coupled  with  RCAS  and  WOPWOP- 
PSU  has  been  developed.  The  UH60A  calculations  are  in  good  agreement  with  experimental 
data.  The  high-order  spectral  difference  code  demonstrates  its  potential  for  capturing  the  rotor 
wake.  Additional  validation  of  codes  for  dynamic  stall  and  blade-vortex  interaction  are  still 
required. 
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